
numADCSamples = 256; % number of ADC samples per chirp
numADCBits = 16; % number of ADC bits per sample
numRX = 4; % number of receivers
numLanes = 2; % do not change. number of lanes is always 2
isReal = 0; % set to 1 if real only data, 0 if complex data0
n_chirps=128;%每一帧中chirps数，no of chirp loops 
frame=1;%第多少帧
numFrames = 5; %总帧数

c=3.0e8;  
B=768e6;       %调频带宽
K=29.982e12;       %调频斜率
T=B/K;         %调频周期
Tc=160e-6;     %chirp总周期
fs=1e7;       %采样率
f0=77e9;       %初始频率
lambda=c/f0;   %雷达信号波长
d=lambda/2;    %天线阵列间距

N = 256;       %距离向FFT点数
M = 128;       %多普勒向FFT点数
Q = 180;       %角度FFT
%% read file
% read .bin file
fileName='D:\20241022\4.bin';
fid = fopen(fileName,'r');
adcData = fread(fid, 'int16');
% if 12 or 14 bits ADC per sample compensate for sign extension
if numADCBits ~= 16
    l_max = 2^(numADCBits-1)-1;
    adcData(adcData > l_max) = adcData(adcData > l_max) - 2^numADCBits;
end
fclose(fid);
fileSize = size(adcData, 1);
% real data reshape, filesize = numADCSamples*numChirps
if isReal
    numChirps = fileSize/numADCSamples/numRX;
    LVDS = zeros(1, fileSize);
    %create column for each chirp
    LVDS = reshape(adcData, numADCSamples*numRX, numChirps);
    %each row is data from one chirp
    LVDS = LVDS.';
else
    % for complex data
    % filesize = 2 * numADCSamples*numChirps
    numChirps = fileSize/2/numADCSamples/numRX;%总chirps数
    numframe=numChirps/n_chirps;%总帧数 no of frames，总chirps数除以每帧chirps数。
    LVDS = zeros(1, fileSize/2);
    %combine real and imaginary part into complex data
    %read in file: 2I is followed by 2Q
    counter = 1;
    for i=1:4:fileSize-1
        LVDS(1,counter) = adcData(i) + sqrt(-1)*adcData(i+2); 
        LVDS(1,counter+1) = adcData(i+1)+sqrt(-1)*adcData(i+3); 
        counter = counter + 2;
    end
        % create column for each chirp
        LVDS = reshape(LVDS, numADCSamples*numRX, numChirps);
        %each row is data from one chirp
        LVDS = LVDS.';
end

    adcData = zeros(numRX, numChirps * numADCSamples);  
    for row = 1:numRX  
        for i = 1:numChirps  
            adcData(row, (i - 1) * numADCSamples + 1:i * numADCSamples) = LVDS(i, (row - 1) * numADCSamples + 1:row * numADCSamples);  
        end  
    end  
  
    % 返回接收器数据  
    adcData; % 接收的总数据 

data_radar_1 = reshape(adcData(1,:),numADCSamples,numChirps);   %RX1
data_radar_2 = reshape(adcData(2,:),numADCSamples,numChirps);   %RX2
data_radar_3 = reshape(adcData(3,:),numADCSamples,numChirps);   %RX3
data_radar_4 = reshape(adcData(4,:),numADCSamples,numChirps);   %RX4
% data_radar=[];            
data_radar(:,:,1)=data_radar_1;     %三维雷达回波数据
data_radar(:,:,2)=data_radar_2;
data_radar(:,:,3)=data_radar_3;
data_radar(:,:,4)=data_radar_4;

%% 3维FFT处理

%距离FFT
range_win = hamming(numADCSamples);   %加海明窗
doppler_win = hamming(numChirps);
range_profile = [];
for k=1:numRX
   for m=1:numChirps
      temp=data_radar(:,m,k).*range_win;    %加窗函数
      temp_fft_1=fft(temp,N);    %对每个chirp做N点FFT
      range_profile(:,m,k)=temp_fft_1;
   end
end

%多普勒FFT
speed_profile = [];
for k=1:numRX
    for n=1:N
      temp=range_profile(n,:,k).*(doppler_win)';    
      temp_fft_2=fftshift(fft(temp,M));    %对rangeFFT结果进行M点FFT
      speed_profile(n,:,k)=temp_fft_2;  
    end
end

%角度FFT
angle_profile = [];
for n=1:N   
    for m=1:M   
      temp=speed_profile(n,m,:);    
      temp_fft_3=fftshift(fft(temp,Q));    %对2DFFT结果进行Q点FFT
      angle_profile(n,m,:)=temp_fft_3;  
    end
end

%% 绘制2维FFT处理三维视图
figure;
speed_profile_temp = reshape(speed_profile(:,:,1),N,M);   
speed_profile_Temp = speed_profile_temp';
[X,Y]=meshgrid((0:N-1)*fs*c/N/2/K,(-M/2:M/2-1)*lambda/Tc/M/2);
mesh(X,Y,(abs(speed_profile_Temp))); 
xlabel('距离(m)');
ylabel('速度(m/s)');
zlabel('信号幅值');
title('2维FFT处理三维视图');
xlim([0 (N-1)*fs*c/N/2/K]);
ylim([(-M/2)*lambda/Tc/M/2 (M/2-1)*lambda/Tc/M/2]);

%% 计算峰值位置
% angle_profile=abs(angle_profile);
% peak=max(angle_profile(:));
% [row,col,pag]=ind2sub(size(angle_profile),find(angle_profile==peak));

% 设置峰值阈值  
peakThreshold = 10;  
 value = db(abs(temp_fft_1)./max(abs(temp_fft_1)))+15;
 a=max(value);
% 使用 findpeaks 函数找到峰值及其位置  
[peakValues, peakLocations] = findpeaks(value,"MinPeakHeight",peakThreshold,'SortStr','descend');  
  
% 初始化目标数组来存储满足条件的峰值位置  
targets = peakLocations;  
  
% 输出检测到的目标（尖峰）位置  
for i = 1:length(targets)  
    disp(['目标', num2str(i), '位置：', num2str(targets(i)/2),'m', ', 值：', num2str(peakValues(i)-15),'db']);  
end  
  
% 输出检测到的目标数量  
disp(['总共检测到的尖峰数量为：', num2str(length(targets)), '个']);

%% 计算目标距离、速度、角度
fb = ((row-1)*fs)/N;        %差拍频率
fd = (col-M/2-1)/(M*Tc);    %多普勒频率
fw = (pag-Q/2-1)/Q;         %空间频率
R = c*(fb-fd)/2/K;          %距离公式
v = lambda*fd/2;            %速度公式
theta = asin(fw*lambda/d);  %角度公式
angle = theta*180/pi;
fprintf('目标距离： %f m\n',R);
fprintf('目标速度： %f m/s\n',v);
fprintf('目标角度： %f°\n',angle);

% %% 绘制目标三维位置信息（雷达为原点）
% % 已知参数  
% phi = pi/2;   % 仰角（单位：弧度）  
%   
% % 将球坐标转换为笛卡尔坐标  
% x = R * sin(phi) * cos(angle);  
% y = R * sin(phi) * sin(angle);  
% z = R * cos(phi);  
%   
% % 创建三维坐标系  
% figure;  
% hold on;  
% axis equal;  
% grid on;  
% xlabel('X');  
% ylabel('Y');  
% zlabel('Z');  
% title('目标位置');  
%   
% % 绘制原点  
% plot3(0, 0, 0, 'ro', 'MarkerSize', 10, 'MarkerFaceColor', 'r');  
% text(0, 0, 0, 'O', 'VerticalAlignment', 'bottom', 'HorizontalAlignment', 'right');  
%   
% % 绘制目标位置  
% plot3(x, y, z, 'bo', 'MarkerSize', 10, 'MarkerFaceColor', 'b');  
% text(x, y, z, 'T', 'VerticalAlignment', 'bottom', 'HorizontalAlignment', 'left');  
%   
% hold off;


